notebook.community

Edit and run



In [ ]:

    
# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import h5py
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
pl.style.use('apw-notebook')
%matplotlib inline
from scipy import interpolate
from scipy.misc import logsumexp

from astroML.utils import log_multivariate_gaussian



In [ ]:

    
DM = 15.6



In [ ]:

    
iso_filename = "/Users/adrian/projects/globber/data/ngc5897/dartmouth_iso_ps1.dat"
iso = ascii.read(iso_filename, header_start=8)



In [ ]:

    
pl.figure(figsize=(4,8))

pl.scatter(iso['gP1']-iso['iP1'], iso['iP1']+DM, alpha=0.5)

pl.xlim(0,1.)
pl.ylim(21,17.5)



In [ ]:

    
interp = interpolate.interp1d(iso['iP1']+DM, iso['gP1']-iso['iP1'], kind='cubic')



In [ ]:

    
itrp_i = np.linspace(17., 21.5, 128)
itrp_gi = interp(itrp_i)



In [ ]:

    
pl.figure(figsize=(4,8))

pl.scatter(itrp_gi, itrp_i, alpha=0.5)

pl.xlim(0,1.)
pl.ylim(21,17.5)



In [ ]:

    
x = np.arange(0,1,0.02)
y = np.arange(17.5,21.,0.04)
shp = (y.size, x.size)
xygrid = np.vstack(list(map(np.ravel,np.meshgrid(x, y)))).T



In [ ]:

    
xymodel = np.vstack((itrp_gi,itrp_i)).T
xymodel.shape



In [ ]:

    
h = 0.03
V = np.diag([h]*2)**2
# W = np.array([[1, -1],   # g-i
#               [0, 1]])  # i magnitude

# # each covariance C = WCW^T
# V = np.einsum('mj,jk->mk', W, V)
# V = np.einsum('lk,mk->ml', W, V)
V



In [ ]:

    
ll = log_multivariate_gaussian(xygrid[:,np.newaxis], xymodel[np.newaxis], V)
ll = logsumexp(ll, axis=1)
ll.shape



In [ ]:

    
pl.figure(figsize=(4,8))

pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
              np.exp(ll).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)

Check that contrast is highest near turnoff



In [ ]:

    
XCov_filename = "/Users/adrian/projects/globber/data/ngc5897/XCov_med.h5"



In [ ]:

    
with h5py.File(XCov_filename, "r") as f:
    bg_X = f['control']['X'][:]
    bg_Cov = f['control']['Cov'][:]
    
    idx = (bg_X[:,0] >= 17.) & (bg_X[:,0] <= 21.5) & (bg_X[:,2] >= -0.5) & (bg_X[:,2] <= 1.5)
    bg_X = bg_X[idx]
    bg_Cov = bg_Cov[idx]
    
    bg_X = bg_X[::50,[2,0]]
    bg_Cov = bg_Cov[::50,[2,0]][:,:,[2,0]]



In [ ]:

    
pl.figure(figsize=(4,8))

pl.plot(bg_X[:,0], bg_X[:,1], alpha=0.4, marker='.', ls='none')

pl.xlim(0,1.)
pl.ylim(21,17.5)



In [ ]:

    
bg_h = 0.025
bg_V = np.diag([bg_h]*2)**2
W = np.array([[1, -1],   # g-i
              [0, 1]])  # i magnitude

# each covariance C = WCW^T
bg_V = np.einsum('mj,jk->mk', W, bg_V)
bg_V = np.einsum('lk,mk->ml', W, bg_V)
bg_V



In [ ]:

    
_V = bg_Cov + bg_V[np.newaxis]
bg_ll = log_multivariate_gaussian(bg_X[np.newaxis], xygrid[:,np.newaxis], _V[np.newaxis])
bg_ll = logsumexp(bg_ll, axis=1)
bg_ll.shape



In [ ]:

    
pl.figure(figsize=(4,8))

pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
              np.exp(bg_ll).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)

The comparison!



In [ ]:

    
pl.figure(figsize=(4,8))

pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
              np.exp(ll - bg_ll).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)



In [ ]: